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Abstract 

Large-scale molecular dynamics simulations are performed to predict the structural 
and thermodynamic properties of liquid krypton using a potential energy function 
based on the two-body potential of Aziz and Slaman plus the triple-dipole Axilrod- 
Teller (AT) potential. By varying the strength of the AT potential we study the 
influence of three-body contribution beyond the triple-dipole dispersion. It is seen 
that the AT potential gives an overall good description of liquid Kr, though other 
contributions such as higher order three-body dispersion and exchange terms cannot 
be ignored. 

The knowledge of interactions in noble gases remains a fundamental question that is not 
completely solved. Despite the simplicity of their closed-shell electronic structure, it is well- 
known that a simple pair potential, though giving the essential features of the structural and 
thermodynamic properties, is not sufficient for a quantitative description, and many-body 
effects have to be taken into accountS. Significant advances have been made when it has been 
demonstratedSil that the static structure factor S{k) at small wave-number, fc, is directly 
related to the long range part of the effective potential between pairs of atoms. It was 
therefore recognized that precise measurements of S{k) could provide a direct observation 
of the details of the interactions. During the past few years, high precision experimentsSH^, 
performed in the range 0.5 < /c < 4 nm~^ using small angle neutron scattering facilities, have 
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confirmed the presence of an additional interaction at large distance that can be associated 
to a three-body contribution at least for Kii and XeQ. 

Long range interactions in noble gases arising from induced dipoles are known as the 
London dispersion forcesi, and are expected to behave as r~^, where r is the interatomic 
distance. Quantum electromagnetic effects acting at very large separationjisl, called retar- 
dation effects, are described by the Casimir-Polder potential^, which is expected to fall off 
as and to have a negligibly small influencei on S{k). It is not surprising that a simple 
Lennard- Jones (12-6) potential is a satisfactory effective pair potential at first sight since 
it can be seen as including mean effects coming from the different dipole-dipole, multipole- 
multipole as well as higher order terms. However, a better understanding of the interactions 
lies in a careful examination of the various genuine contributions of the multipole expansion, 
which contains two-body dd, dq, qq terms, etc., as well as irreducible three-body contribu- 
tions such as ddd, ddq, dqq, qqq, etc., where d and q denote the dipole and quadrupole 
moment, respectively. According to Barker and Henderson§, four -and more- body terms 
are very small and can be neglected. The two-body potential of noble gases is often taken in 
the very accurate Hartree-Fock-dispersion form propounded by Aziz and Slaman0'El, while 
three-body dispersion effects are often represented by the ddd term in the form given by 
Axirod and Tellei0 that represents the major contribution. The other three-body contri- 
butions mentioned above, which can be modelled by the expressions of Belli, are known to 
have a small influence although it is not clear whether the agreement with the experiment 
could be improved when these are taken into account. At high densities, beside dispersion 
terms, three-body overlap contributions such as exchange effects may come into playM. 

It is tempting to consider an effective three-body potential of the Axirod- Teller^ (AT) 
form 

1 + 3 cos 9i cos 9j cos 9^ 

U^iVi, Tj, n) = V -3-3-3 (1) 

ij ik jk 

to take all these three-body contributions into account by modifying the strength u to 
include them empirically. In Eq. ([l|) 9i, 9j and 9k denote, respectively, the angles at vertex 
i, j and k of the triangle k) with sides Vij = \rj — rj|, rj^ = \rk — rj| and rjk = jr^ — rj\. 
Usually, the AT potential is designed to represent the ddd contribution only with a value 
u = Uddd = 2.204 X 10~^^ J nm^ as prescribed by Leonard and Barker^. Such an effective 
three-body potential was recently used to study the liquid-vapor phase equilibria of argonill 
as well as the small-A; part of the static structure factor of krypton in the dense liquidi. For 
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instance, to interpret correctly their experimental data, Guarini et alE have increased the 
strength u by 65 %, indicating that other three-body contributions play an important role. 
Conversely, in a recent workil, we have shown that the missing contributions would have the 
effect of slightly reducing the strength of the AT potential. Nevertheless, these studies are 
not completely conclusive since calculations where carried out within the integral equations 
theory, in which the three-body contribution has to be treated as an approximate state- 
dependent effective pair potentialil, which is not unique. 

The purpose of this article is to analyze the influence of the three-body interactions taken 
in the form of Eq. (|l]), to decide what is the effective value of the strength u for an accurate 
description of liquid Kr. Three different situations are considered by varying u: a value of 
z/ = that corresponds to the pure two-body potential given by Aziz and Slamanlll, a value 
of z/ = Uddd representing only the ddd contribution, and a value Ueff = ^-Q^^ddd which is 
believed to reproduce the critical parameters of Kr within 1 %i. Changing the strength of 
the AT potential is a convenient means to measure of the departure from the ddd term and, 
by comparison with the experiments, could provide useful information on the contribution 
coming from all other three-body terms that are omitted in the modelling of the interaction 
between Kr atoms. This work represents an extension to the liquid phase of our preceding 
studies on Kr in the gaseous phaseSil. Thus, for a given interaction model, we carry out 
molecular dynamics (MD) simulations to determine S{k) as well as the internal energy and 
the virial pressure. In this context, MD is a powerful methodiH since the three-body 
interactions can be tackled without incurring the shortcomings of the approximate integral 
equations. As a matter of fact, the treatment of the three-body potential (|l|) is not subject to 
arbitrariness since the resulting forces are derived in an exact manne , as for the two-body 
ones. 

In order to extract a meaningful structure factor S{k) from the MD simulations, the 
pair correlation g{r) is the key quantity that has to be calculated as precisely as possible. 
Therefore, we have performed large-scale MD in the sense that (i) a large enough simulation 
cell has been considered so that g{r) has a sufficient spatial extension to yield a correct S{k) 
by Fourier transform, especially at low k, and (ii) a large number of time steps are produced 
in order to get a significant part of the phase space trajectory, essential for the statistics. 
As three-body forces are involved, this represents a huge amount of computer time, and we 
have used a parallel algorithm described in some details in a previous work@. Simulations 
with the different values of u in the AT potential have been performed in the microcanonical 
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(NVE) ensemble with a time step At = 5.67 • 10~^^ s using = 6912 particles in a cubic 
cell subject to the standard boundary conditions. The cutoff radius of the interactions is 
Tc = 2.5rm, where = 0.4008 nm is the minimum of the AS potential. For the three-body 
potential, this implies that triplets of particles in which two or three distances of separations 
are greater than Tc are ignored in the calculation of the forces. In order to investigate the 
influence of Vc for the two-body potential, a value of Ar^ is taken in specific cases. The typical 
duration of the runs is 113 ps from which 1300 independent configurations are extracted for 
the statistical analysis of the physical quantities. Six states of liquid Kr from the vicinity of 
the critical point to that of the triple point have been studied, which are those investigated 
by Guarini et a/.i and by Barocchi et al^, respectively at low- and large-A;. These states 
correspond to temperature T = 199 K and densities n = 12.10, 11.66 and 11.31 nm~^, 
T = 169 K with n = 14.57 and 14.22 nm'^, and T = 130 K with n = 16.83 nm'^. 

In Fig. 1, we present the large k behavior of S{k) for the three different temperatures 
and we compare the MD curves, calculated with the three different values of u, to the 
experimental data of Barocchi et a/.0. As the temperature increases and density decreases, 
the first sharp diffraction peak as well as the subsequent oscillations become less pronounced. 
This is well predicted by the MD results since a remarkable agreement with the experiments 
is found whatever the values of u. It can be seen that the three-body contributions have 
only a minor influence even if a more careful examination shows that, at low temperature, 
the height of the first peak is better described without the AT potential while, at high 
temperature, a closer agreement is obtained when it is included. At low-k, the three potential 
energy functions give rise to completely different behaviors of S{k) and the best results seem 
to be those obtained with the ddd interaction. 

Let us now focus on the low-A; part of S{k) in more details. We compare the MD 
curves to the recent small-angle scattering data of Guarini et alE along the T = 199 K and 
T = 169 K isotherms in Fig. 2 (a) and (b), respectively. At the highest temperature, near 
the critical isotherm, a good agreement is found with the experiment at the three densities 
with the ddd potential. Examining the influence of u at density n = 12.10 nm~^, it appears 
that the AS two-body potential alone {u = 0) is not able to predict the small-/c part of 
S{k) satisfactorily, while the effective AT contribution (z/ = i^eff) gives rise to a structure 
factor which is underestimated. The chain curve corresponds to S{k) obtained with the AS 
potential alone and a cutoff radius = Avm- The influence of is seen only below 2.5 nm~^ 
and the curvature is slightly changed. As a result, this would have the effect of increasing the 
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values of S{k) in this region whatever the potential used. Taking rc = 4rm when the three- 
body contribution is considered is computationally too costly and could represent a challenge 
for future developments of the present MD code. However, for the three-body potential, it 
is worth noting that, even with Vc = 2.5rm, configurations in which one pair of particles of 
a triplet is separated up to 5r.m are taken into account in the calculation of the forces. In 
addition, the AT potential given by Eq. (|I|) decays as therefore taking a cutoff radius 
larger than 2.5rm seems not to be necessary. The same observation can be drawn along the 
T = 169 K isotherm, as shown in Fig 2 (b), and in this case a very good agreement with 
the experimental curves is seen when the ddd potential is taken into account. Interesting 
enough, both experimental data sets of Guarini et a/.i and Barocchi et a/.0, which connect 
to each other very well around 4 nm~^, lie between the curves corresponding to u = Vddd 
and V = Veffi for k values in the range between 2.5 and 6.5 nm~^. 

At T = 130 K the situation is less clear as it can be seen in Fig. 3. While the AS 
potential is also not sufficient to predict the small-fc behavior of S{k), this time the best 
concordance with the experiment is obtained by using v = Veff- Nevertheless, the latter 
effective AT contribution would have a tendency to underestimate the PVT data while the 
ddd strength gives the best prediction. In this case, taking again a cutoff radius of Arm for 
the two-body potential alone (chain curve) has only a negligible influence on S{k), therefore 
our results for the two- plus three-body might be correct. It should be stressed that the 
amplitude of S{k) is very small at such a low temperature and high density state, and 
the relative difference between the MD curves with Vddd and z/e// is of the same order of 
magnitude than the dispersion of the experimental data points, which is about 10 %. At 
this stage, we refrain from drawing any conclusion and it would be desirable to dispose of 
an accurate small angle scattering experiment for this thermodynamic state. 

We also examine the influence of the three-body potential on the internal energy E and 
the virial pressure P gathered in Table 1. Long-range corrections have been applied, due to 
the truncation of the AS and the AT potentials during the simulation at r^. = 2.5rm. The 
corrections to the two-body part of the energy and the pressure are respectively — 0.153?T,r^ 
and —0.312 {nr'^)^. For the three-body parts, the corrections are estimated numerically 
using the integral equation methodil. For the energy, we obtain 0.007 (nrf^)'^ for Uddd and 
0.012 (^r^)^ for i^eff, while for the pressure we get 0.031 (nrf^)^ for Uddd and 0.051 (nr^)^ 
for i^eff- These quantities are expressed in the unit of the minimum of the AS potential 
e = 2.777 ■ 10^^^ J. The influence of u on the energy is moderate and does not exceed 9 % 
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with respect to that of the two-body potential, even with J^eff- On the contrary, it has a large 
effect on the pressure. For the two-body potential alone, the pressure is always negative, 
while by including the three-body potential contribution, it becomes positive in the majority 
of cases. In addition, z/g// yields pressure values which are too high, while the ddd strength 
gives the best predictions, even if it is always smaller compared to the experimental data. 
Again, the calculated pressures with Uddd and z/g// enclose the experimental values. 

Regarding the results presented above, it appears that the details of the interaction model 
has no significant influence on S{k) at large k, as shown in Fig. 1. Even the AS two-body 
potential alone is able to reproduce the structure factor with a good degree of accuracy. For 
the structure factor at small scattering angle and the pressure, the three-body interactions 
cannot be ignored in the liquid state. Moreover, it is seen that the triple-dipole contribution 
gives the best agreement with the experiments and therefore represents the main three-body 
effect. A more precise examination shows that the structural and thermodynamic properties 
depart substantially from the experiments even with the model combining the AS two-body 
potential plus ddd contribution. Indeed, discrepancies remain (i) on the small k part of 
S{k) in the range between 2.5 and 6.5 nm~^, where the calculated values are higher that 
the experimental ones, especially for isotherms T = 169 K and T = 130 K, and (ii) on the 
pressure where the theoretical results are systematically below the measurements, and even 
remain negative for two thermodynamic states. Since the pressure varies linearly with z/, as 
it can be seen in Table 1, better results should be obtained by increasing the value of u^ff 
between 1.20i'iidd and 1.25i'ddd, whatever the temperature. 

We are led to the same conclusion as Guarini et a/.i that the value of u in Eq. (P has 
to be increased with respect to that of the triple-dipole Uddd- However, by the light of the 
present MD calculations, an effective strength u^ff = l.Qbuddd seems to be too important, 
and we estimate that the additional three-body terms beyond the ddd one should represent 
up to 25 % of it. Now the question arises to know what is the nature of the missing 
contributions that will take a non negligible part in the interaction model. According to 
Copeland and KestneJll who studied liquid argon, two majors three-body potentials beyond 
the ddd one play an important role, namely the exchange overlap and the ddq dispersion 
acting respectively at short and long interatomic distances. Both potentials are known to 
have significant infiuenceS'i, however, while the ddq potential has the same sign as the 
ddd contribution, the exchange one has an opposite sign. Therefore, it will be of primary 
importance to investigate their interplay in liquid krypton and whether they improve the 
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description of the interactions in hquid Kr when added to the ddd term. It has been recently 
shown by van der Hoef and Maddenll in their simulation study of liquid argon that the ddq 
contribution on the pressure is small but not negligible and it would be interesting to extend 
these results in the case of Kr, not only for the pressure but also for the structure factor. 

At very small wave-number, i.e. k < 2.5 nm~^, the MD results of S{k) underestimate the 
experimental curves as well as the PVT data, which is particularly visible along the isotherm 
T = 199 K displayed in Fig. 2 (a). As MD simulations are concerned, we are unavoidably 
confronted to finite size effects and this fact might be attributed mainly to two factors: the 
truncation of the pair-correlation functions at the half of the box size, and the use of a cutoff 
radius of the interactions. Therefore, the S{k) calculated by Fourier transform are subject 
to large uncertainties and must be taken cautiously. Moreover, near the critical region, the 
correlation length can exceed the size of the simulation box and for this reason the present 
MD simulations might be not able to catch the correct behavior of S{k) at T = 199 K. As 
a matter of fact, as pointed out by Wilding^ and Roverell, in the thermodynamic limit, 
critical phenomena, like the divergence of S'(O), are smeared out and shifted. 

In conclusion, MD simulations have been carried out for Kr in the liquid phase for which 
new small angle scattering experiments were recently performediS. This study completes 
preceding works on the low density and high temperature states of where it was 

demonstrated that the Aziz and Slaman two-body potential associated to the Axirod and 
Teller triple-dipole contribution gives a excellent representation of interactions in the gaseous 
phase. Here, we have shown that the latter potential energy function predicts the essential 
characteristics of structural properties in liquid Kr, even though we believe that an accurate 
description of S{k) at low k and the thermodynamic properties requires that other three- 
body contributions such as the dipole-dipole-quadrupole and the exchange overlap potentials 
are taken into account. Large scale MD simulations including these additional contributions 
will be performed in the near future. 
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Captions 



Figure 1. Structure factor S{q) for T = 130 K, at n = 16.83 nm-^, for T = 169 K, at 
n = 14.57 nm~^ and for T = 199 K, at = 12.10 nm~^ from the top to the bottom 
(the curves for T = 169 K and T = 130 K are shifted upwards by an amount of 
1 and 2, respectively), calculated by molecular dynamics with u = (dashed hues) 
= t^ddd (sohd hues) and u = Ueff (dotted hues) as described in the text. Open circles 
correspond to the experimental data of Ref. ^ while full circles stand for the PVT 
data of Ref. 

Figure 2. Structure factor S{q) at small scattering angle for isotherm (a) T = 199 K, 
at n = 12.10 nm~^, n = 11.66 nm~^ and n = 11.31 nm~^, and (b) T = 169 K, at 
n = 14.57 nm~'^.and 14.22 nm~^. The curves for n = 11.66 nm~'^ and n = 11.31 nm~^ 
are shifted upwards by an amount of 0.5 and 1, respectively, and that of 14.22 nm~^ 
by an amount of 0.2. Molecular dynamics results are carried out with u = (dashed 
lines), u = Vddd (solid lines) and v = v^ff (dotted lines) as described in the text. The 
chain curve corresponds to MD results with u = and a cutoff radius of 4rm. Open 



circles correspond to the experimental data of Ref. 30, open triangles correspond to 



the experimental data of Ref. ^ and full circles stand for the PVT data of Ref. |33| . 

Figure 3. Structure factor S{q) at small scattering angle for T = 130 K, at = 16.83 
nm^^. Same captions as in Fig.l. 

Table 1. Excess internal energy, E'^^/Ne, and pressure, P (r^)^ /e, calculated by molecular 
dynamics for the different thermodynamic states considered in this work. Subscripts 
2 and 3 stand respectively for the two- and two- plus three-body parts of the internal 
energy. Pexp corresponds to the experimental values of Ref. ^ The numbers in brackets 
represent the standard deviations that affect the last decimal of the temperature, 
energy and pressure extracted from the simulation. 
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Table 



T{K)p 


(nm ^) 


T^/T^ddd 


Ef/Ne 


Ef/Ne 


Prl/e 


P (bar) Pexp 


(bar) 


130 (1) 


16.83 





-4.962 (6) 


— 


-0.87 (4) 


-375.4 








1 


-4.953 (6) 


-4.689 (6) 


-0.12 (4) 


-51.8 









1.65 


-4.958 (6) 


-4.522 (6) 


0.33 (4) 


142.4 




169 (1) 


14.57 





-4.151 (7) 


— 


-0.37 (5) 


-159.6 








1 


-4.127 (7) 


-3.938 (7) 


0.06 (4) 


25.9 


61.7 






1.65 


-4.126 (7) 


-3.818 (7) 


0.30 (4) 


129.4 




169 (1) 


14.22 


1 


-4.034 (7) 


-3.855 (7) 


-0.05 (4) 


-21.5 


20.1 


199 (1) 


12.10 





-3.392 (8) 


— 


-0.06 (4) 


-25.9 








1 


-3.358 (9) 


-3.232 (9) 


0.11 (4) 


49.6 


73.1 






1.65 


-3.344 (8) 


-3.138 (8) 


0.27 (4) 


116.5 




199 (1) 


11.66 


1 


-3.245 (8) 


-3.127 (8) 


0.09 (4) 


38.8 


55.7 


199 (1) 


11.31 


1 


-3.165 (8) 


-3.053 (8) 


0.06 (4) 


25.9 


46.3 



Table 1 
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k (nm'^) 




k (nm"^) 




k (nm'^) 



